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Abstract 

Be^"'"(aq) hydrolysis can to lead to the formation of multi-beryllium clusters, but 
the thermodynamics of this process has not been resolved theoretically. We study 
the hydration state of an isolated Be^"^ ion using both the quasi-chemical theory of 
solutions and ab initio molecular dynamics. These studies confirm that Be^^(aq) is 
tetra-hydrated. The quasi-chemical approach is then applied to then the deproto- 
nation of Be(H20)4^+ to give BeOH(H20)3+. The calculated pK^ of 3.8 is in good 
agreement with the experimentally suggested value around 3.5. The calculated en- 
ergetics for the formation of [Be • OH • Be]^^ are then obtained in fair agreement 
with experiments. 

Key words: beryllium, ab initio molecular dynamics, deprotonation, speciation, 
quasi-chemical theory 



1 Introduction 

Beryllium metal has properties that make it technologically very attractive[l], 
but these advantages are severely counterbalanced by the high toxicity of 
inhaled beryllium dust, which causes chronic beryllium disease in a subset of 
exposed individuals [2] . The etiology of this autoimmune disease [3] is poorly 
understood, but the final disease state is characterized by lung failure. 

Aqueous beryllium chemistry is also incompletely understood. Experiments 
suggest that beryllium mediated hydrolysis of water leads to the formation of 
multi-beryllium clusters [1,4]. This same mechanism is likely involved in the 
dissolution of Be^^(aq) and of importance in environmental clean-up strate- 
gies. Such mechanisms likely underlie deposition of beryllium in biomaterials, 

* Corresponding author: Phone: 505-667-8624; Fax: 505-665-3909 
Email address: lrp@lanl.gov (Lawrence R. Pratt). 



Preprint submitted to Elsevier Science 



2 February 2008 




Fig. 1. Structures representing the deprotonated tetra-aquo cation, BeOH[H20]3"'', 
and the aggregate formed by the coalescence of one such unit and a Be[H20]4^"'", 
with explusion of a water molecule. 

and thus perhaps in the development of chronic beryllium disease. Thus a 
molecular understanding of Be^"'"(aq) and formation of multi-beryllium species 
would provide a foundation for addressing these issues of wide importance. 
This letter takes an initial theoretical step in understanding the aggrega- 
tion/disaggregation of beryllium clusters in water. 

Early electronic structure calculations on beryllium hydration were performed 
on small clusters [5,6] and some attempted to include the second hydration 
shell in terms of the reaction field approach [7] or explicitly [8]. Molecular 
dynamics simulation [9] showed that assuming pair-wise intermolecular in- 
teractions lead a hydration number was six (6) whereas including three-body 
effects brought the hydration number down to four (4) consistent with solution 
X-ray diffraction experiments. 

Be^"*" hydration has also been considered within the Car-Parrinello approach 
[10,11,12]. There a hexa-hydrate cluster quickly dissociated to give a tetra- 
hydrate structure [10]. The ab initio molecular dynamics approach for a Be^"*" 
atom in a box of 31 water molecules also showed that the ion was tetra- 
hydrated [12]. That work also suggested an effect of the second hydration 
shell water molecules on some of the bond lengths for the central Be(H20)4^^ 
structure, noted further below. 

Those earlier works have addressed important issues about Be^+ solvation, 
but the thermodynamical aspect of solvation has not been considered specif- 
ically. We address the solvation thermodynamics by employing two distinct 
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theoretical approaches. We study the solvation structure of Be^"^ using ab 
inifAo molecular dynamics methods and confirm a stable tetra-coordination. 
The quasi-chemical theory of solutions then provides a further analysis why 
the tetra-hydrate is the most stable species. This approach has been used be- 
fore to address Fe^''"(aq) speciation [13], and the hydration of H+ [14], H0~ 
[15], Li"*" [16], and Na"*" [17] ions. As mentioned above, beryllium solution 
chemistry holds the particular challenge of the formation of multi-beryllium 
clusters. The present studies on Be(H20)4^"'" deprotonation and on formation 
of [Be ■ OH ■ Be]^"*" provide first steps in describing that aggregation process on 
a molecular basis, as suggested by Fig. 1. 



2 Ab initio moleculcir dynamics 

The ab initio molecular dynamics (AIMD) simulations were performed with 
the VASP program [18,19]. The simulation system comprises 32 water molecules 
and one Be^"*" ion. The box length was set to 9.71 A based on the experimental 
partial molar volume of Be in water [20] . Vanderbilt ultrasoft pseudopotentials 
[21,22] were used to describe the core-valence interactions for all the atoms. 
The valence orbitals were expanded in plane waves with a kinetic energy cut- 
off of 29.1 Ry. All the hydrogen atoms were replaced by deuterium and an 
integration timestep of 0.5 fs was chosen. 

The initial configuration was obtained by placing a Be^+ ion in a bath of 
water molecules. That configuration was energy minimized before initiating 
the AIMD simulations. At the end of the classical energy minimization, the 
Be^"*" ion was found to be pcnta-hydratcd (based on R(BeO) < 2.5 A). In an 
exploratory AIMD run (data not reported), the initial configuration for the 
AIMD simulation was obtained from a classical molecular dynamics simula- 
tion. In that case the ion was hexa-hydrated, but it too quickly reverted to 
the tetra-hydrated form. 

In the first 1.0 ps of the AIMD simulation the temperature was maintained 
at 300 K by scaling the velocities. After this initial phase, a microcanonical 
(NVE) ensemble simulation was performed for about 2.5 ps. 

Within the first 160 fs of the ca. 3.5 ps of the AIMD simulation, the coordina- 
tion number changed to four (4) and stayed so for the rest of the simulation. 
The energy in the NVE simulation was -472.25±0.08 (2cr) eV, suggesting good 
energy conservation. The mean temperature was 316.2±21.6 K. 

Fig. 2 shows the oxygen(water) radial distribution around the beryllium ion. 
Observe that the inner shell is physically sharply defined and the hydration 
number is four (4). Further structural characterization of Bc[Il20]4^"'" solvation 
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Fig. 2. Radial distribution of oxygens around Be^^. The right ordinate (dashed hue) 
gives the coordination number, n(r). 

is collected in table 1. Earlier[12] it was suggested that the second hydration 
Table 1 

Solvation structure of Be[H20]4^^. Inner: Results for the inner-shell hydration; 
Outer: Results for the outer-shell hydration. Cluster: Results pertaining to a gas- 
phase cluster (or isolated molecule) obtained with the methods in section III. All 
lengths are in A and angles are in degrees. Values in parentheses are for the isolated 
II2O molecule. 



AIMD 





Inner 


Outer 


Cluster 


R(BeO) 


1.64±0.06 




1.64 


ZOBeO 


109.3±5.9 




111.4±5.9 


ROH 


1.02±0.04 


1.00±0.03 


0.98±0.01 (0.97) 


ZHOH 


108.5±5.6 


106.5±5.5 


111.4 (105.7) 



shell has a significant influence on the OH bond length of the inner shell water 
molecule. In the present case, we do find a slight lengthening of the OH bond 
in the inner-shell water, but the statistical uncertainties here and in the earlier 
study [12], suggest that the difference between the inner-shell and outer-shell 
water is subtle. The increase in the HOH angle for the inner shell water is 
similar to the earlier study [12] and is in line with the values obtained for an 
isolated Be(H20)4^+ cluster. Also notice that the oxygen atoms in the inner- 
shell are nearly tetrahedrally distributed around the central Be^"*" ion, but that 
structure does fluctuate somewhat. 
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3 Quasi-chemical Theory 



In the quasi-chemical theory [23], the region around the solute of interest is 
partitioned into inner and outer shell domains. In the present study, the inner 
shell, where chemical effects are important, is treated quantum mechanically. 
The outer shell contributions have been assessed using a dielectric continuum 
model. In principle, a variational check of this partition is available (see [17]). 

Hydration state of Be^+: The inner shell reactions pertinent to Be^"^ hy- 
dration are: 



Be2+ + nHaO # Be(H20)n^+ 

The free energy change for these reactions were calculated using the Gaussian 
programs [24]. The Be • (H20)n^''' (n = 0. . .6) clusters were geometry optimized 
in the gas phase using the B3LYP hybrid density functional [2 5] and the 6- 
31-|-G(d,p) basis set. Frequency calculations confirmed a true minimum, and 
the zero point energies were computed at the same level of theory. Single point 
energies were calculated using the 6-311+G(2d,p) basis set. A purely inner- 
shell n = 5 cluster could not be obtained; the optimization gave structures 
with four (4) inner and one (1) outer sphere water molecule. For n — 6 both a 
purely inner-shell configuration, and a structure with four (4) inner and two 
(2) outer shell water molecules were obtained. The quasi-chemical theory here 
utilizes only the inner-shell structure. 

For estimating the outer shell electrostatic contribution, the ChelpG method 
[26] was used to obtain partial atomic charges. Then with the radii set de- 
veloped by Stefanovich et al.[27], surface tessera were generated [28], and the 
solvation free energies of the clusters were calculated using a dielectric contin- 
uum model [29]. With this information and the binding free energies for the 
chemical reactions, a primitive quasi-chemical approximation to the excess 
chemical potential of Be^''"(aq) in water is: 



/5/^Be2+(a,)~-ln(i^nPH20" 

where = KJ^'^ exp -/? (AiBe(H20)„2+ - Wrzo) ■ -^n^°^ is the equihbrium 
constant for the reaction in an ideal gas state, with n of Eq. 1 the hydra- 
tion number of the most stable inner shell cluster, and P — l/keT. The den- 
sity factor PH2O appearing in eq. 1 reflects the actual density of liquid wa- 
ter and its effect is accounted for by including a replacement contribution of 
— nkeT ln(1354). A detailed statement on standard states and this replacement 
contribution can be found in Grabowski et al. [14]. 
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Fig. 3. Quasi-chemical contributions of the solvation free energy of Be^^(aq). Specif- 
ically, the plot (left ordinate) is -k^T x In (^KnpU20^^ vs. n predicted by the 
primitive quasi-chemical approximation; see Ref. [16]. n is the number of inner 
shell water molecules surrounding the anion. A: —RThiKn^ — ni?rin(1354); x: 
^Be(H20) 2+ ^ "'/^H20- An observation volume of radius 2.0 A centered on the cation 
defined the inner shell. Using a smaller radii did not make an appreciable difference. 

From figure 3 it is clear that the tetra-aquo cation is the most stable form in 
solution. This is also consistent with the predictions of the AIMD simulations. 
As fig. 3 indicates, neglecting solvation effects would have forced us to conclude 
that both the tetra- and hexa-hydrates should be observed. However, it is the 
substantial unfavorable solvation of the hexa-hydrate that precludes its pres- 
ence in the solvent. The calculated solvation free energy of -567.7 kcal/mole is 
in the right range and is in good agreement with the value (-574.6 kcal/mole) 
cited in [20] . Solvation free energy values for these highly charged species may 
have significant uncertainties, and the agreement here should not mask the real 
difficulties facing ion-solvation thermodynamics. The agreement in inner-shell 
hydration number between AIMD simulations and quasi-chemical calculations 
is, however, a non-trivial result. 

In table 2 some of the energies relevant to the n — A case above is collected in 
addition to other free energy values that we will use in analyzing the pKa of 
the tetra-aquo cluster. 

pKa of Be(H20)4^+: The acidity of Be(H20)4^+ is described by the 



[BeOH(H20)3+] [H^ 



[Be(H20)4 



2+1 



(2) 



corresponding to the reaction 



Be(H20)4'+ ^ BeOH(H20)3+ + H+ (3) 
under standard conditions (1 M ideal solution). The reaction as written re- 
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Table 2 

Electronic energy (a.u.), thermal corrections (a.u.) to the free energy, and ex- 
cess chemical potential (kcal/mole) using dielectric continuum approximation with 
charges obtained at B3LYP/6-311+G(2d,p). 





E 










-U.UloUO 




Be(H20)42+ 


-320.12296 


0.07259 


-230.9 


BcOH(H20)3+ 


-319.92375 


0.05970 


-73.9 


Be20H(H20)63+ 


-563.42051 


0.13225 


-412.0 


Be20H=^+ 


-103.83704 


-0.00627 




H2O 


-76.45951 


0.00298 


-7.7 


HO- 


-75.82779 


-0.00771 




HO • (H20)3- 


-305.32036 


0.04705 


-72.3 



quires us to know the hydration free energy of the proton, estimates of which 
have large uncertainties [14]. Alternatively, using a reference reaction such as 
the dissociation of water can obviate the need for knowing the proton hydra- 
tion free energy. Thus consider 

Be(H20)4^+ + HO- ^ BeOH(H20)3+ + H2O (4) 

Knowing the energetics of eq. 4 permits calculation of the pKa of the tetra- 
aquo cation, according to 

Ka = Kx{Kj[il,0]) . (5) 

Here K is the equilibrium ratio of cq. 4 with all concentration units the same, 
Ky, = [H0~] [H+] the standard ion product for water; here pK^ = 15.7 [30]. 
An added advantage of using a reference reaction as above is that some can- 
cellation of errors can be encouraged. 

The excess chemical potential of HO~(aq) in reaction 4 above is also ob- 
tained within the quasi-chemical approach and with the tri-hydrated quasi- 
component [15]. The computed value of -105 kcal/mole is in good agreement 
with recently reported values [31]. 

If all standard concentrations are 1 M, then using the values in tabic 2 the 
free energy change for reaction 4 is -18.5 kcal/mole. Eq. 5 then yields pKa ~ 
3.8. 
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By fitting experimental free energy clianges for tlie case of low total Be^"*" 
concentration [1,4], it is found that Be(H20)4^+ exists in appreciable amounts 
only below a pH of 3.5. The present calculated pKg, is in excellant agreement 
with these observations. This value of pKa has the standard interpretation 
that the deprotonated complex BeOH(H20)3'*" is above a thousand times more 
probable than Be(H20)4^''" at neutral pH. That a spontaneous deprotonation 
is not observed in our simulations is reflective of the limited simulation time 
and the possibility of this process being activated. Nevertheless, when an 0H~ 
was introduced into the AIMD simulation by extraction of a proton distant 
from the beryllium ion, ligand exchange by proton shuffling was accomplished 
in less than a ps. 

Formation of [Be • OH • Be]^+: The [Be ■ OH • Be]^+ cluster is one of the 
many clusters that beryllium forms [1], but it is the simplest. Other clus- 
ters could be constructed with this unit. Thus understanding the formation 
thermodynamics of this cluster is of first interest. 

The solvation free energy of the complex [Be • OH ■ Be]^"*" is obtained from the 
following reaction. 



Here is is assumed that each of the originally tetra-hydrated Be^+ loses one 
water and gains a H0~ to form the complex (fig. 1). Using eq. 1 for the present 
case, the quasi-chemical estimate of the solvation free energy of [Be-OH-Be]^'*' 
is obtained as -834.5 kcal/mole. 

With this solvation free energy estimate, we can enquire about the thermody- 
namics of the complexation reaction below. 



The change in the excess chemical potential for the reaction is 405.9 kcal/mole, 
whereas the ideal contribution i.e., for an the ideal gas at 1 atm pressure, 
is -424.2 kcal/mole for this change. This gives the net free energy change 
of -18.3 kcal/mole not accounting for standard concentration. Converting to 
the standard concentration of 1 M adjusts this by — 2-RTln24.46 to the 
calculated free energy change of -22.1 kcal/mole. 

An experimental value for the complexation reaction is -14.4 kcal/mole [1,4]. 
The calculated energetics are in the right range, but the agreement is only fair. 

Note also that a small difference between large numbers is being computed; 
thus even minor errors will tend to get amplified. A physical conclusion is that 
the hydration contribution is more 20 times larger than the net standard free 
energy change for this reaction; neglecting hydration effects would lead to a 



Be • OH ■ Be^+ + 6H2O # Be • OH ■ Be ■ (H20)6 



3+ 



(6) 



2Be^+ + HO- # [Be • OH • Be]=^+ 



(7) 
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qualitatively incorrect result. Below we consider ways to improve upon these 
initial estimates. 



4 Concluding Discussions 

The quasi-chemical approach leads to free energies that are in reasonable 
agreement with available experimental estimates. But this agreement should 
not obscure the severe approximations that have been made in applying the 
theory to practical calculations. For divalent cations, it is unclear if purely a 
inner-shell complex would provide adequate accuracy of solvation free energy 
estimate. There are potentially two possible avenues for improvement which 
we indicate below. 

First, the outer-sphere hydration contributions can be obtained using classical 
molecular mechanics approaches instead of the dielectric continuum model 
used above. Such an approach is now being undertaken for water clusters on 
monovalent cations (Asthagiri, et al. in preparation) and this can likely be 
used for Be(H20)4^+ as well. 

Second, it is possible to expand the chemical potential of the tetra-aqua cation 
in terms of its own inner- shell (i.e. the Be2+ cation's 2"^ shell). This seems par- 
ticularly natural here because the inner shell structures are physically definite 
(fig. 2). The hydration contribution of those shell structures could be ob- 
tained using a dielectric continuum model. Also it is possible to apply a lower 
level of quantum chemical approximation to describe the coupling of the I''* 
and 2"^^ shells. Then with an improved estimate for the solvation free energy 
of the tetra-aqua cation, one could obtain a better estimate for the solva- 
tion free energy of Be^+. This approach is the iterated quasi-chemical scheme 
and a variant has been successfully applied to hard-sphere fluids [32] . In the 
present case, however, this approach would be daunting for many reasons. (A) 
Obtaining statistically representative 2"°'-shell complexes in the gas-phase is 
difficult. (B) The rigid- rotor harmonic oscillator approximation is dubious for 
applications to such clusters. (C) Applying the quasi-chemical theory to such 
large clusters demands consideration of packing effects. This latter issue is 
actively under study at present (Ashbaugh and Pratt, in preparation). 

The neglect of packing aspect is certainly one the reason underlying the 
merely fair agreement for the [Be • OH ■ Be]^+ thermochemistry. It is also likely 
that anharmonic effects are non-negligible for clusters such as Be20H^+ and 
Be20H(H20)6^'^. But despite these substantial limitations, it is heartening to 
obtain the qualitatively correct trends. In our pKa calculations, however, by 
arranging the equation to have similarly sized species on both sides of the 
equality we were able to mitigate these uncertainties. 
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